1/1 


AD-A157  736  ITERATIVE  AND  FADE  SOLUTIONS  FOR  THE  WATER-WAVE 

DISPERSION  RELATIONS)  COASTAL  ENGINEERING  RESEARCH 
CENTER  VICKSBURG  MS  H  S  CHEN  ET  AL.  APR  85  CERC-85-4 

F/G  20/4 


UNCLASSIFIED 


NL 


AD-A157  796 


ttw 


fl 


US  Army  Corps 
of  Engineers 


[■•r-v 


MISCELLANEOUS  PAPER  CERC-85-4 

ITERATIVE  AND  PADE  SOLUTIONS  FOR  THE 
WATER-WAVE  DISPERSION  RELATION 

by 

H.  S.  Chen  and  E.  F.  Thompson 
Coastal  Engineering  Research  Center 

DEPARTMENT  OF  THE  ARMY 
Waterways  Experiment  Station,  Corps  of  Engineers 

PO  Box  631 

Vicksburg,  Mississippi  39180-0631 


April  1985 
Final  Report 

Approved  For  Public  Release;  Distribution  Unlimited 


s 


DTIC 

ELECTE| 
JUL1  9  866  ■ 


Prepared  for  DEPARTMENT  OF  THE  ARMY 
US  Army  Corps  of  Engineers 
Washington,  DC  20314-1000 

under  Waves  at  Entrances  Work  Unit  31673 


85  6  28  122 


& 


m 


Destroy  this  report  when  no  longer  needed.  Do  not  return 
it  to  the  originator. 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated 
by  other  authorized  documents. 


I  itcmlon  Tor  / 
rifflS  ORAM  bT 

•  dtic  tab  Tj 


GRAAI 
DTIC  TAB 
Unannounced 
Justification- 


Distribution/ 


Availability  Codes 


(Avail  and/or 


The  contents  of  this  report  are  not  to  be  used  for 
advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an 
r  ficial  endorsement  or  approval  of  the  use  of 
such  commercial  products. 


_ Unclassified _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 

REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER 


Miscellaneous  Paper  CERC-85-4 


4.  TITLE  (end  Subtitle) 


ITERATIVE  AND  PADE  SOLUTIONS  FOR  THE  WATER-WAVE 
DISPERSION  RELATION 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

ECIPIENT'S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  A  PERIOD  COVERED 


Final  report 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHOR!*; 

H.  S.  Chen 
E.  F.  Thompson 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

US  Army  Engineer  Waterways  Experiment  Station 

Coastal  Engineering  Research  Center 

PO  Box  631,  Vicksburg,  Mississippi  39180-0631 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

DEPARTMENT  OF  THE  ARMY 
US  Army  Corps  of  Engineers 

Washington,  DC  20314-1000 _ 


B.  CONTRACT  OR  GRANT  NUMBER!*; 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Waves  at  Entrances  Work 
Unit  31673 
___________  _ 

April  1985 

13.  NUMBER  OF  PAGES 
18 


14.  MONITORING  AGENCY  NAME  B  ADDRESS!//  dllletent  Iron,  Controlling  Olllce)  15.  SECURITY  CLASS,  (of  thle  report) 

Unclassified 

IS*.  DECLASSIFICATION/DOWNGRADING 
SCHEDULE 


l«.  DISTRIBUTION  STATEMENT  (a!  thle  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (ot  the  ebetrecl  entered  In  Bloc*  20,  II  dlllerent  Iron,  Report) 


IB.  SUPPLEMENTARY  NOTES 

Available  from  National  Technical  Information  Service,  5285  Port  Royal  Road, 
Springfield,  Virginia  22161. 

19.  KEY  WORDS  (Ccnttnue  on  rererme  aide  #/  neceeemry  and  Identity  by  block  number) 

-Water  waves — Dynamics 
Dispersion  relations-  ( LGy- 
Pade'  approximant^  / 

Iterative  methods  (Mathematics)  A  — 

20k  ABSTRACT  CCeerttmue  on  to  roe  ee  .M  R  nmeeeeeer  end  Identity  *7  Block  nienber) 

'  — - y-tfhe  iterative  and  Pade^olutions  of  the  propagating  and  evanescent  wave 

modes  of  the  dispersion  relation  for  the  linear  water-wave  theory,  are.  given  in 
this  report.  -Their  corresponding  FORTRAN  subroutines  are  also  provided.  The 
Pade^ solution  is  generally  more  efficient  than  the  iterative  one  in  most  engi¬ 
neering  applications,  /it.  f  *  rl  s  -  — - 


DO  i'ZTt.  1473  EDITION  OF  t  MOV  SB  IS  OBSOLETE 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dele  Entered) 


yf 


PREFACE 


This  report  is  a  by-product  of  the  Waves  at  Entrances  Work  Unit ,  Harbor 
Entrances  and  Coastal  Channels  Program,  Civil  Works  Research  and  Development, 
at  the  Coastal  Engineering  Research  Center  (CERC)  of  the  US  Army  Engineer 
Waterways  Experiment  Station  (WES). 

The  report  was  prepared  by  Drs.  H.  S.  Chen  and  E.  F.  Thompson,  Coastal 
Oceanography  Branch  (COB),  CERC,  under  direct  supervision  of  Dr.  J.  R.  Houston 
Chief,  Research  Division,  CERC,  and  under  general  supervision  of  Dr.  R.  W. 
Whalin,  Chief,  CERC. 

The  authors  wish  to  acknowledge  the  assistance  of  Ms.  Willie  A.  Brown, 
Physicist,  CERC,  in  performing  some  of  the  calculations. 

Commander  and  Director  of  WES  during  the  publication  of  this  report  was 
COL  Robert  C.  Lee,  CE.  Technical  Director  was  Mr.  F.  R.  Brown. 


CONTENTS 


PREFACE  .  1 

PART  I:  INTRODUCTION .  3 

PART  II:  APPROXIMATE  SOLUTIONS  OF  THE  DISPERSION  RELATION  .  4 

Propagating  Wave  Mode .  4 

Evanescent  Wave  Mode .  7 

PART  III:  CONCLUSION .  15 

REFERENCES .  16 

APPENDIX  A:  LISTING  OF  SUBROUTINES  .  A1 


2 


ITERATIVE  AND  PAPE  SOLUTIONS  FOR  THE 
WATER-WAVE  DISPERSION  RELATION 


PART  I:  INTRODUCTION 

1.  In  linear  water-wave  theory,  the  dispersion  relation 

2 

u)  =  gk  tanh  kh  (1) 

describes  the  relationship  between  wave  frequency  u>  =  2n/T  and  wave  number 
k  =  2rt/L  for  gravity  waves  on  water  depth  h  where  g  is  the  gravitational 
acceleration,  T  the  wave  period,  and  L  the  wavelength.  It  has  been  shown 
(Mei  1983)  that  for  a  given  u)  and  h  Equation  1  has  a  pair  of  real  roots, 
k  =  ±kQ  ,  representing  the  propagating  wave  mode  and  an  infinite  number  of 
discrete  imaginary  roots,  k  =  ±i<*Q  (n  =  1,2 . »)  ,  representing  the  eva¬ 

nescent  wave  mode.  The  terms  k  and  a  are  real  and  i  =  J-l  .  The  prop- 

o  n 

agating  wave  mode  represents  progressive  waves,  while  the  evanescent  wave  mode 
represents  local  disturbances  or  standing  waves  which  decay  out  in  the  far 
field.  The  latter  is  particularly  important  in  the  study  of  wave- structure 
interaction  which  involves  measurement  or  calculation  in  the  near  field,  such 
as  wave  forces  on  a  structure,  flow  field  and  free  surface  elevation  near  a 
structure,  wave  maker  design,  etc.  Since  there  is  no  exact  explicit  solution 
of  Equation  1  for  k  ,  the  wave  number  is  usually  obtained  by  an  iterative 
method.  Hunt  (1979)  derived  an  approximate  solution  for  the  propagating  wave 
mode  kQ  by  the  Pade  approximation.  Hunt's  solution  was  given  also  in  a 
Coastal  Engineering  Technical  Note  (CETN)  (US  Army  Engineer  Waterways  Experi¬ 
ment  Station,  Coastal  Engineering  Research  Center  1982)  which  is  currently 
being  revised.  The  approximate  solution  is  very  accurate,  even  in  the  short 
wave  range.  In  this  study  a  Pade'  solution  for  the  evanescent  wave  mode 
is  also  developed.  Additionally,  comparisons  of  accuracy  and  central  pro¬ 
cessing  unit  (CPU)  time  for  the  iterative  and  Pade'  solutions  are  presented 
as  well  as  FORTRAN  subroutines  for  calculating  wave  number  by  both  schemes. 


PART  II:  APPROXIMATE  SOLUTIONS  OF  THE  DISPERSION  RELATION 


2.  Equation  1  is  rewritten  in  the  following  form: 


y  =  kh  = 


a 

tanh  kh 


(2) 


2 

where  a  =  ui  h/g  .  Then  a  pair  of  real  roots  ±kQ  exists  as  solutions  of 
Equation  1.  The  solutions  can  be  represented  graphically  as  the  intersection 
points  of  the  two  curves,  y  =  kh  and  y  =  o/tanh  kh  ,  as  shown  in  Figure  1. 


using  Equation  2  and  following  the  dashed  arrow  lines  to  the  intersection 
point  as  illustrated  in  Figure  1.  A  FORTRAN  subroutine,  PKITER,  for  calcu¬ 
lating  kQ  using  the  iterative  approximation  is  provided  in  Appendix  A. 
Pade  approximation 

4.  By  employing  the  Pade  scheme,  Hunt  (1979)  obtained  an  approximate 
solution  for  Equation  1  in  the  form  of  the  following  rational  function: 


5.  The  first  nine  values  of  are:  =  0.66667  ,  d =  0.35550  , 

d3  =  0.16084  ,  d4  •*  0.06320  ,  d5  =  0.02174  ,  dfe  =  0.00654  ,  d?  =  0.00171  , 
dg  =  0.00039  ,  dg  =  0.00011  .  A  FORTRAN  subroutine,  PKPADE,  for  calculating 
using  the  Pade  approximation  is  provided  in  Appendix  A. 

6.  Calculated  results  of  k  for  various  values  of  ct  by  the  itera- 

°  -8  -6 

tive  scheme  with  the  truncation  errors  of  10  and  10  and  the  Pade  scheme 
with  M  =  6  (six  terms)  and  M  =  9  (nine  terms)  are  given  in  Table  1.  It 

Table  1 

Comparison  of  the  Results  of  the  Iterative  and  Pade 
Approximations  for  the  Propagating  Wave  Modes 


Calculation 

Technique 

PKPA6* 

PKPA9* 

PKIT6** 

PKIT8** 


0.001  3. 16280480E-02 

3. 16280480E-02 
3.16275484E-02 
3.16280430E-02 


Relative 

Error 

1.58E-07 

1.58E-07 

-1.56E-05 

0.00 


PKPA6 

PKPA9 

PKIT6 

PKIT8 


0.010  1 . 00166971E-01 

1 . 00166971E-01 
1.00167469E-01 
1 . 00166968E-01 


2.99E-08 

2.99E-08 

5.00E-06 

0.00 


_ (Continued) _ _ _ _ 

*  PKPA6  and  PKPA9  are  the  Pade  solutious  of  M  =  6  and  M  =  9,  respectively . 

**  PKIT6  and  PKIT8  are  the  iterative  solutions  with  truncation  errors  of  10 
-8 

and  10  ,  respectively.  The  values  of  PKIT8  are  used  as  reference  values 

for  calculating  the  relative  errors. 
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Table  1  (Concluded) 


Calculation 

lr 

Relative 

Technique 

CT 

O 

Error 

PKPA6 

0.050 

2.25487360E-01 

-3.10E-08 

PKPA9 

2 . 25487360E-01 

-3.10E-08 

PKIT6 

2 . 25487841E-01 

2.10E-06 

PKIT8 

2 . 25487367E-01 

0.00 

PKPA6 

0.100 

3 . 21595936E-01 

8.40E-08 

PKPA9 

3 . 21595936E-01 

8.40E-08 

PKIT6 

3 . 21596366E-01 

1.42E-06 

PKIT8 

3 . 21595909E-01 

0.00 

PKPA6 

0.500 

7 . 71706799E-01 

5.81E-06 

PKPA9 

7 . 71704463E-01 

2.78E-06 

PKIT6 

7 . 71702707E-01 

5.07E-07 

PKIT8 

7.71702316E-01 

0.00 

PKPA6 

1.000 

1 . 19985791E+00 

1.49E-04 

PKPA9 

1 . 19968005E+00 

1.18E-06 

PKIT6 

1 . 19967884E+00 

1.67E-07 

PKIT8 

1 . 19967864E+00 

0.00 

PKPA6 

2.000 

2 . 06858341E+00 

1.57E-03 

PKPA9 

2.06522702E+00 

-5.38E-05 

PKIT6 

2.06533819E+00 

2.42E-08 

PKIT8 

2.06533814E+00 

0.00 

PKPA6 

3.000 

3.0I938365E+00 

1.63E-03 

PKPA9 

3.01459484E+00 

3.72D-05 

PKIT6 

3 . 01448277E+00 

-3.32E-09 

PKIT8 

3.01 448278E+00 

0.00 

PKPA6 

3.500 

3.51050861E+00 

1.20E-03 

PKPA9 

3.50657600E+00 

7.62E-05 

PKIT6 

3.50630882E+00 

0.00 

PKIT8 

3 . 50630882E+00 

0.00 

PKPA6 

4.000 

4.00588653E+00 

8.04E-04 

PKPA9 

4 . 00298773E+00 

7.93E-05 

PKIT6 

4.00267030E+00 

0.00 

PKIT8 

4 . 00267030D+00 

0.00 

PKPA6 

5.000 

5 . 00205761E+00 

3.21E-04 

PKPA9 

5 . 00067224E+00 

4.37E-05 

PKIT6 

5 . 00045361E+00 

0.00 

PKIT8 

5 . 00045361E+00 

0.00 

PKPA6 

6.000 

6 . 00082 146E+00 

1.25E-04 

PKPA9 

6.00017538+00 

1.69E-05 

PKIT6 

6 . 00007372E+00 

0.00 

PKIT8 

6.00007372E+00 

0.00 

6 


indicates  that  the  Pade  solutions  of  M  =  6  and  M  =  9  are  within  accuracies 

of  0.1  percent  and  0.01  percent,  respectively,  in  comparison  with  the  itera- 

“8 

tive  solutions  with  truncation  error  of  10  .  Besides  the  convenience  of  an 

explicit  form,  the  Pade  solution  in  Equation  3  also  saves  CPU  time.  For  cal¬ 
culating  the  twelve  values  of  a  ,  the  amount  of  CPU  time  saved  is  shown  below: 

PKPA6*  PKPA9  PKIT6  PKIT8 

less  than 

0.0012  0.0012  0.9302  0.9998 


*  Time  in  seconds. 

Evanescent  Wave  Mode 

7.  In  Equation  1  there  are  imaginary  roots  k  =  ia  where  a  is  real. 
Inserting  this  expression  for  k  into  Equation  1  leads  to  the  following 
equation: 

2 

ui  =  -g a  tan  ah  (4) 


8.  The  real  roots  of  Equation  4  correspond  to  the  imaginary  roots  of 
Equation  1.  Equation  4  can  be  rewritten  as 


y  = 


P 

ah 


=  tan  ah 


(5) 


2 

where  a  =  (<u  h)/g  as  before. 

9.  These  roots  are  the  intersection  points  of  the  curves  y  =  -a/ah  and 
y  =  tan  ah  .  Since  tan  ah  has  an  infinite  number  of  branches,  there  is  an 
infinite  number  of  discrete  roots,  a  =  (n  =  1,2,...,°°)  ,  as  shown 

graphically  in  Figure  2. 


Iterative  approximation 


10.  An  iterative  algorithm  for  calculating  an  is  established  by  using 
Equation  5  and  following  the  dashed  arrow  lines  to  the  intersection  point  as 
illustrated  in  Figure  2.  A  FORTRAN  subroutine,  EKITER,  for  calculating  a^ 


using  the  iterative  approximation  is  provided  in  Appendix  A.  Calculated  re- 

-8  -4 

suits  of  a  (n  =  1,2,..., 10)  with  truncation  errors  of  10  and  10  are 
n 

given  in  Tables  2  and  3,  respectively. 


I  W-y=-§r 

Figure  2.  Imaginary  roots  of  the  dispersion  relation 


Pade  approximation 

11.  For  brevity,  Equation  U  can  be  rewritten  as 

o  =  -x  tan  x  (6 

where  x  =  orh  .  In  order  to  take  appropriate  account  of  the  multiplicity  of 
roots,  let 

x  =  nn  -  x’  ,  -  -  <  x'<  0  and  n  =  1,2,...,°°  (7 

Substituting  Equation  7  into  Equation  6  gives 

ct  =  (nn  -  x')  tan  x'  (8 

From  Equation  8,  we  have 


x'  -*  —  as  ct  -*■  0 

nn 


.  n 

x'  -*■  -  as  ct  -»  « 


>2 


12.  Therefore,  for  convenience  the  following  rational  function  for 
approximation  is  chosen: 


(10) 


13.  Note  that  Equation  10  is  by  no  means  the  only  form  for  the  approxi¬ 
mate  solution  of  Equation  8,  as  one  can  readily  understand  from  the  Pade  ap¬ 
proximation  (Baker  1975  and  Bender  and  Orszag  1978).  Expanding  the  functions 
of  Equation  8  and  Equation  10  into  a  series  and  matching  the  coefficients  as 

in  the  Pade  approximation,  we  obtain  the  following  values  of  b  (m  =  1,2,3, 

m 

4,5): 


bl 

b2 

b3 

b4 

b5 

n  =  1 

0.20264 

0.06159 

0.01395 

0.00182 

-0.00065 

n  =  2 

0.10132 

0.01283 

0.00083 

-0.00003 

-0.00001 

n  =  3 

0.06755 

0.00532 

0.00017 

-0.00001 

0.0 

n  =  4 

0.05066 

0.00289 

0.00006 

0.0 

0.0 

n  =  5 

0.04053 

0.00181 

0.00003 

0.0 

0.0 

n  =  6 

0.03377 

0.00124 

0.00001 

0.0 

0.0 

n  =  7 

0.02895 

0.00090 

0.00001 

0.0 

0.0 

n  =  8 

0.02533 

0.00068 

0.00001 

0.0 

0.0 

n  =  9 

0.02252 

0.00054 

0.0 

0.0 

0.0 

n  =  10 

0.02026 

0.00043 

0.0 

0.0 

0.0 

and  the  solution  of  the  evanescent  wave  number  is 


(ID 


14.  We  should  note  that  computational  experiments  in  this  study  indi¬ 
cate  that  an  increase  in  M  (such  as  M  =  6,  7,  8,  9,  or  10)  does  not  improve 
accuracy  in  the  short  wave  range  at  a  =  4  ,  5  ,  and  6  .  Therefore,  we  choose 
M  =  5  in  the  calculations.  A  FORTRAN  subroutine,  EKPADE,  for  calculating 
an  for  n  up  to  10  using  Equation  11  is  provided  in  Appendix  A.  Calculated 


results  of  a  (n  =  1,2,..., 10)  are  given  in  Table  4. 
n 

15.  Comparisons  of  accuracy  of  the  results  in  Tables  2,  3,  and  4  are 

shown  in  Table  5.  It  indicates  that  the  Pade  solution  in  Equation  11  gives  an 

approximation  with  maximum  error  of  1.5  percent  at  a  =  6  in  comparison  with 

-8 

the  iterative  solution  with  a  truncation  error  of  10  .  The  error  becomes 

smaller  as  ct  becomes  smaller.  We  should  note  that  the  large  value  of  a 
corresponds  to  short  waves;  for  example,  T  is  about  3  ~  1  sec  for 
a  =  4  ~  6  ,  at  water  depth  h  =  10  m  .  In  this  short  wave  range,  most 
engineering  applications  are  relatively  insignificant.  The  CPU  time  savings 
is  about  threefold  or  sixfold  that  of  the  iterative  solutions,  depending  on 
the  truncation  errors  accepted.  The  CPU  time  (in  seconds)  for  calculating 
the  ten  modes  (M  =  1,2,..., 10)  of  the  twelve  values  of  a  is  shown  below: 

EKIT  4*  EKIT  8*  EKPA** 

0.03662  0.07324  0.01221 

*  EKIT8  and  EKIT4  stand  for  the  iterativg 
solution^  with  truncation  errors  of  10" 
and  10"  ,  respectively.  The  values  of 

EKIT8  are  used  as  reference  values  for 
calculating  the  relative  errors. 

**  EKPA  stands  for  the  solution  of  Equation 

11. 

16.  All  calculations  were  conducted  on  a  Control  Data  Cyber  170  Model 
760  computer.  Without  loss  of  generality,  the  water  depth  h  was  taken  to 
be  1.0  in  the  calculations. 


Table  4 

✓ 

Pade  Solution  (of  Equation  11)  with  M  =  5  for  the  Evanescent  Wave  Mode 
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Table  5 


Comparisons  of  Accuracy  of 

the  Iterative 

and 

Pade  Approximations 

for 

the  Evanescent 

Wave 

Mode 

Relative  Error 

CT 

EKIT8 

EKPA 

0.001 

0 

0 

3.2E-09 

0.010 

1.3E-07 

0 

6.7E-08 

0.050 

3.2E-09 

0 

1.3E-08 

0.100 

1.3E-08 

0 

1.5E-07 

0.500 

4.8E-07 

0 

6.4E-06 

1.000 

2.0E-06 

0 

7 . 8E-05 

2.000 

1.4E-05 

0 

-5.1E-04 

3.000 

7.9E-06 

0 

-5.7E-03 

3.500 

7.3E-06 

0 

-8.8E-03 

4.000 

5.6E-06 

0 

-1.1E-02 

5.000 

1.5E-05 

0 

-7.0E-03 

6.000 

7.7E-06 

0 

1.5E-02 
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PART  III:  CONCLUSION 


17.  The  Pade  solutions  of  the  dispersion  relation  in  Equation  1,  given 
by  Equations  3  and  11,  are  rational  functions  which  are  preferable  for  handy 
calculation.  Their  accuracies--being  within  0.1  percent  for  Equation  3  and 
1.5  percent  for  Equation  11 — and  CPU  time  savings--being  tremendous  for  Equa¬ 
tion  3  and  at  least  threefold  for  Equation  11 — should  enhance  wave¬ 
scattering  calculations  in  most  engineering  applications.  Four  FORTRAN  sub¬ 
routines  calculating  wave  numbers  for  the  propagating  and  evanescent  wave  modes 
by  iterative  and  Pade  schemes  are  documented  in  Appendix  A. 
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APPENDIX  A:  LISTING  OF  SUBROUTINES 


1*  In  this  appendix,  the  following  FORTRAN  subroutines  for  calculating 
wave  number  are  documented: 

a.  PKITER:  calculating  the  propagating  wave  mode  by  the  iterative 
approximation. 

b.  PKPADE:  calculating  the  propagating  wave  mode  by  the  Pade 
approximation  in  Equation  3  with  M  =  9  . 

c.  EKITER:  calculating  the  evanescent  wave  mode  by  the  iterative 
approximation. 

d.  EKPADE:  calculating  the  evanescent  wave  mode,  n  =  1  to  10,  by 
the  Pade  approximation  in  Equation  11  with  M  =  5  . 

2 

2.  In  the  subroutines,  the  inputs  are  as  follows:  OMGG  is  u  /g  ;  H  is 

til 

water  depth  h  ;  TOUR  is  truncation  error;  N  is  the  n  evanescent  wave 

mode;  and  PI  =  n  -  3.141592654.  The  outputs  are  PK  ,  the  propagating  wave 

number  k  ,  and  EK  ,  the  evanescent  wave  number  a 
o’’  n 

3.  The  listing  of  the  FORTRAN  programs  is  as  follows: 


SUBROUTINE  PK  I  TER  ( OMGG ,  H * PK ,  T0L.R ) 


0MGG=0MGA**2/G*  H=WATER  DEPTH *  PK«UAVE  NO, 


Y=0MGG*H 
XJ=Y 
10  XI =XJ 

XJ=Y/TANH(XI ) 

IF( ABS(XI-XJ) .GT.TOLR)  GOTO  10 

PKH=XJ 

PK=PKH/H 

RETURN 

END 

SUBROUTINE  PKPADE ( OMGG , H  *  PK ) 


0MGG=0MGA**2/G*  H=WA TER  DEPTH*  PK=WAVE  NO. 


Y=0MGG*H 

XX= Y* ( Y+ 1 . / ( 1 . +Y* ( 0 . 66667+ Y* < 0 . 35550+ Y* ( 0 , 16084+ Y* ( 0 , 06320+Y* 

1  ( 0 . 02174+Y* ( 0 . 00654+  Y* ( 0 . 001 71+ Y*  < 0 , 00039+ Y*0 .00011))))*));}) 

X=SQRT(XX) 

PK-X/H 

RETURN 

END 

SUBROUTINE  EK I TER  < OMGG , H  >  EK  *  N  *  P I  *  T0LR ) 

0MGG=0MGA**2/G*  H=WATER  DEPTH*  EK=EUANESCENT  WAVE  NO.*  N^NTH  MODE. 


Y=0MGG*H 
XJ=FL0AT(N)fPI 
DX=XJ 
10  XI^XJ 

XJ-ATAN(-Y/XI)+DX 

IF(AFS/XI-XJ).GT.TOLR»XJ)  GOTO  10 

EKH=XJ 

EK-EKH/H 

RETURN 

END 


A1 


SUBROU TINE  EKPABE ( OMGG *H*EK*N*PI) 

C  .  . . . . - . . 

C  0HGG=0MGA**2/G>  H=WATER  IiEF'TH »  EK-EUAMESCENT  WAVE  NO.*  N-NTH  MODE . 

c - - - 

DIMENSION  8(5*10) 

DATA  B/0. 20264 >0.06159 >0.01395* 0.00182* -0.00065*0. 10132 >0.01283* 
0. 00083*-0, 00003*-0. 00001 *0. 06755*0. 00532*0. 00017*-0. 00001 *0.0* 

0  *  05066  *  0 . 00289 « 0 . 00006  *  0 . 0  *  0 . 0  *  0 ♦ 04053  *0.00181*0. 00003  >0.0  *0.0* 
0.03377*0.00124*0.00001*0.0*0.0*0.02895*0.00090*0.00001*0.0*0.0* 
0.02533*0.00068*0.00001*0.0*0.0*0.02252*0.00054*0.0*0.0*0.0* 
0.02026*0.00043*0.0*0.0*0.0/ 

F'I2=PI/2. 

Y=OMGG*H 
EKH=FLOAT (N)fPI 

EKH=EKH-PI2*(1 .-!,/( l.+YJ|t(B(l  *N)+Y*(B<2*N)+Y#(B(3*N5+Y]!f(B'(  4*  N)  + 

1  Y*B *'  5  *  N )  .»)>))) 

EK=EKH/H 
RETURN 
END 


